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NUMERICAL  SOLUTION  OF  SUPG  FINITE-ELEMENT 
METHOD  FOR  SUPERSONIC  VISCOUS  FLOW 

Xu  Guoqun  and  Zhang  Guofu 

Nanjing  Aeronautical  College,  Nanjing  210016 

Abstract:  A  streamline  upwind/Petrov-Galerkin  (SUPG) 
weighted  residual  formalism  was  developed  in  the  paper  for  the 
quasi-simplified  Navier-Stokes  equations.  Numerical  computations 
were  made  for  Burger's  equation,  the  nonviscous  shock  wave 
reflection  problem,  as  well  as  supersonic  laminar  flow  over  a 
flat  plate  and  compression  corner  flow  by  using  the  method.  The 
results  of  the  computations  show  that  this  method  is  accurate, 
convergent,  and  stable. 

Key  words:  Navier-Stokes  equations,  supersonic  viscous  flow, 
finite-element  method,  SUPG  method. 

I.  Introduction 

When  the  convection  flow  term  is  present,  the  differential 
operator  is  nonself-adjoint;  the  coefficient  matrix  in  the 
Galerkin  method  is  nonsymmetric.  The  method  lacks  the  advantage 
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of  optimal  approximation.  Often  this  situation  causes 
fluctuations  in  the  numerical  solution.  For  problems  involving 
high  Pe  numbers  or  large  Re  numbers,  this  phenomenon  is 
especially  serious.  To  control  convection  of  each  element,  the 
mesh  can  be  made  finer;  however,  in  this  case  a  large  amount  of 
internal  memory  will  be  used.  In  1979,  Hughes  and  Brooks  [l] 
proposed  the  SUP G  method.  The  method  has  better  stability  and 
higher  accuracy,  capable  of  effectively  processing  high-speed 
flow  problems,  including  the  problem  in  which  the  interruption 
plane  lies  in  the  flow  plane.  The  present  authors  verified  that 
the  SUPG  method  and  the  PG  (Penalty  Galerkin)  method  developed  by 
Baker  [2]  are  somewhat  equally  effective  for  high-speed  viscous 
flow  problems.  However,  expansion  and  application  of  the  SUPG 
method  is  more  promising  than  the  PG  method  [3],  In  references 
[4,  5],  the  SUPG  finite-element  method  is  used  to  solve  for  the 
incompressible  Navier-Stokes  equation  set.  In  reference  [6],  the 
Euler  equation  set  is  solved.  However,  up  to  now  the  authors  did 
not  find  any  report  in  which  the  SUPG  finite-element  method  was 
used  in  solving  problems  of  supersonic  viscous  flow. 

II.  Model  Equations 

Below,  Burgers'  equation  is  discussed: 

!r+*£-*ip-0’  -><«<'•»» 

- —  1<X<1 

*(/,-l)-0.5t  w(/,l) - 0.5,  t>  0 

Define  the  weighted  function  in  the  element  as 
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w,  -  N,  +  P,  (2) 

In  the  equation,  is  the  interpolation  function,  a  function, 
which  is  continuous  at  the  element  boundary.  P^  is  the 
perturbation  of  the  weighted  function,  a  function,  which  is 
discontinuous  at  the  element  boundary.  Then,  we  write  out  the 
Petrov-Galerkin  expression  form  of  Eq.  (1)  as 

+ 1  Lp'(!f+ « 

In  the  equation,  ng^  is  the  number  of  elements;  in  this  paper, 
the  linear  interpolation  value  is  taken  as  u.  Therefore,  there 
is  no  second-order  derivative  term  in  the  second  integration  term 
of  Eq.  (3).  In  the  SUPG  finite-element  method,  the  perturbation 
P.  of  the  weighted  function  can  be  taken  as 

P,  —  oA*  Ml  (4) 

ox 

Here,  OC  is  the  algorithm  parameter,  and  At  is  the  duration  of 
the  time  step. 

Use  Eq.  (3)  to  solve  for  the  above-mentioned  problem  of  a 

specific  solution.  Fig.  1  shows  the  results  of  computations  when 
-4 

k  «  10  .  Fig.  1(a)  and  (b)  are,  respectively,  the  results  when 

0  and  when  OC  =  0.3.  The  numerical  solutions  shows 
fluctuations  upstream  and  downstream  of  the  shock  wave;  this  is  a 
nonphysical  solution,  which  fails  to  satisfies  the  entropy- 
increase  condition  [7]  expressed  in  the  second  law  of 
thermodynamics.  Fig.  1(d)  shows  the  solution  when  Of  -  1.0.  The 
solution  appears  as  an  excessive  divergence  in  the  neighborhood 
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of  the  shock  wave,  thus  showing  lower  accuracy.  Fig.  1(c)  shows 
the  results  when  ot=  0.5.  At  this  point,  the  virtual  fluctuation 
of  the  finite-element  solution  in  the  neighborhood  of  the  shock 
altogether  disappears;  in  addition,  accuracy  is  also  ensured. 


<O«-0  (b)  a  —  0.3  (c)a-O.S  (d)  a  - 1.0 


Fig.  1.  Solutions  of  Burgers'  equation 
when  k  =  10-4  given  in  computations 
accurate  solution 

•  —  numerical  solution  as  given  in  the  paper 

Fig.  2(a)  and  (b)  show,  respectively,  the  numerical  results 

-2 

when  k  *  0  and  when  k  =  10  .  From  the  figure,  we  can  see  that 

the  numerical  solution  and  the  accurate  solution  are  in  good 

-4 

agreement.  In  addition,  when  k  is  less  than  or  equal  to  10  , 

the  shock  wave  travels  only  one  mesh  point.  The  method  exhibits 
high  discriminability  for  a  shock  wave. 

III.  Quasi-simplified  Navier-Stokes  Equations  and 
Their  Variation  Form 

Assume  that  t  stands  for  time,  u  and  v  are  components  of  the 
gas  flow  velocity  in  the  x  and  y  axes,  p,  f  ,  T,  and  fi,  are  the 
pressure,  density,  temperature,  and  viscosity  coefficients  of  the 
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gas;  and  v  is  the  adiabatic  index.  To  simplify  the  computations, 
the  second-order  derivative  terms  in  the  x-direction  (along  the 
direction  of  the  object  surface)  can  be  neglected  in  supersonic 
and  hypersonic  problems.  Then  in  the  case  of  dimensionless 
nonsteady  state  for  an  ideal  gas  at  constant  specific  heat,  the 
Navier-Stokes  equation  set  can  be  simplified  as 


dt  dx  dy  dy 


(5) 


Here,  U,  o-  ,  and  are  vectors;  and  A 2  are  the  matrixes 
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Fig.  2.  Calculated  solutions  of  Burgers'  equation 
when  OC  -  0.5 

_ _  accurate  solution  •  -  numerical  solution 

as  given  in  the  paper 


To  bring  the  equation  set  to  a  closure,  the  state  equation 
and  the  viscosity  law  equation  should  be  added.  The  state 
equation  is 


P 


l 

vMl 


PT 


C«> 


There  is  also  the  constant -volume  specific  heat  equation, 

C,  —  i)  .  The  viscosity  equation  adopts  the  Sutherland 


relationship 


--7-VI  1  +  sr/rm 
t  •+■  y/r. 


(7> 


In  the  above-mentioned  equation,  Re*,^  is  the  Reynolds  number  of 
the  incoming  flow,  and  Pr  is  the  Prandtl  number  (taken  as  0.72  in 
the  paper).  When  calculating  ft  ,  S'  ■  114K. 

The  weak  solution  is  taken  for  Eq.  (5).  In  addition, 

Green's  formula  is  used,  thus  we  establish  the  variation  equation 
in  the  following  form: 
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(8) 


UM'fr+A' af+A,fr-+*)_'’  fjM" 

+  s  j  r>  (o  + *■  + A-  fy- + *  )■ " 


N,odx 


In  Eq.  (8),  I  is  the  4x4  element  matrix.  is  the  fundamental 

function  of  the  interpolation  value,  which  is  continuous  at  each 
element  boundary.  In  the  paper,  the  interpolation  values  of  four 
nodal  points  are  taken  for  p,  u,  v,  and  T: 

(9) 


Assume  that  the  natural  boundary  conditions  are  homogeneous 
in  Eq.  (8)  for  spatial  divergence.  We  have  the  semi-divergence 
finite-element  equation 

MU  +  (C  +  F)t/  —  0  (10) 

In  the  equation,  M  *  M(U,t),  which  is  the  mass  matrix. 

C  =  C(U,t),  which  is  the  rigidity  matrix  (or  the  convection 
matrix).  F  is  the  divergence  matrix;  and  U  is  the  derivative  of 
U  with  respect  to  time.  The  matrix  in  Eq.  (10)  is  the 

v 

aggregation  of  each  element. 


IV.  Selection  of  Perturbation  of  Weighted  Function 

In  the  SUPG  finite-element  method,  the  perturbation  P^  of  a 
weighted  function  is  related  to  the  fraction  matrix  of  the  first- 
order  derivative  terms.  Now,  P^  is  written  as 
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In  the  equ  ion,  't  •  A^  ,  and  T  •  A2are,  respectively,  the 
artifi^al  viscosity  coefficients  along  the  x  and  y  directions. 
In  the  paper,  nr  is  taken  as  related  to  the  duration  of  the  time 
step.  Experience  with  computations  proves  that  this  approach  is 
effective. 

r-a.A*  (12) 


In  the  equation,  (X  is  the  algorithm  parameter.  There  are 
different  values  to  be  taken  for  equation  sets  of  different 
kinds.  According  to  reference  [6],  the  steady-state  time-step 
duration  At  of  the  entire  computation  can  be  estimated  by  using 
the  following  equation. 

“■  (not  solving  for  the  sum)  (13) 

In  the  equation,  the  subscript  i  indicates  the  coordinate 
direction;  the  superscript  e  indicates  the  compilation  number 

a 

of  the  elements;  h^.  is  the  length  of  element  e  along  the 
direction  x^  (refer  to  Fig.  3);  a^  is  the  spectrum  radius  of  A 
(the  maximum  absolute  value  for  the  characteristic  value). 
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X  ! 

Fig.  3.  Finite  element  of  equal  parameters 


V.  Estimation  and  Revision  Algorithm 

Write  out  Eq.  (10)  as  the  all-implicit  form 

M'.'haV:;*  +  C Vi .€/.*;{»  -  0  (14) 


In  the  above  equation,  the  subscript  is  the  time  step  and  the 
superscript  is  the  number  of  iterative  substitutions  of  the 
corresponding  time  step  — {/<;♦;>  •  Now  let  us  assume 


A«  -  -  «“i,  (15) 

and  assume  that  f  is  the  relaxation  factor.  Then  Eq.  (14)  can 
be  written  as 


In  the  equation 


M*Aa  -  R 


(16) 


M*  -  M"h  +  +  Fi'i,)  ( 17  ) 

R  -  -  (C'.'i,  +  F«> ,)C/S,  ( 18) 


We  can  see  from  Eq.  (16)  that  if  M*  is  considered  as  mass 
and  R  is  the  driving  force,  then  Ape  corresponds  to  acceleration. 
The  iterative  substitution  of  each  time  step  ensures  the  phase 
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equilibrium  between  inertial  force  and  driving  force.  Also  from 
Eq.  (18)  we  know  that  if  a  convergent  solution  exists  at  a  time 
step,  then  driving  force  R  should  gradually  approach  zero  (the 
residue  approaches  zero)  with  the  process  of  iterative 
substitution;  therefore,  the  acceleration  Aoc  also  gradually 
tends  to  zero. 

By  combining  Eqs.  (15)  and  (16),  we  have  the  following 
estimation-revision  algorithm: 


(l)  i-o  (*  fc&fttttt)  a 


<4)  R  -  -M'lUa.'U  -  Ci'i.f/i ‘It  («£*) 
<5)  M*A a  -  R  ) 


(6)  a.*:: » - 

(7)  U:x?  - 


a'Jl,  +  A  a 


U“l,  +  fiAt&a 


C 


KEY:  a  -  (i  is  the  number  of  iterative  substitution) 
b  -  (estimation  stage) 
c  -  (residual  forced 
d  -  (revision  stage) 

If  a  further  step  of  iterative  substitution  is  required,  the 
computation  returns  to  step  4.  In  this  paper,  the  nonsymmetric 
linear  equation  sets  are  found  by  using  this  algorithm;  the  wave 
matrix  technique  is  used  to  solve  the  equation  sets.  The 
condition  for  convergence  for  each  time  step  is 

max|/?|  <  s,  (19) 

JO-*  ,  Since  there  is  an  estimation  stage,  generally 
speaking  the  revision  times  are  not  many,  only  three  or  four 


times  to  satisfy  Eq.  (19).  A  method  with  too  many  revision  times 
is  not  desirable.  When  the  number  of  revision  times  is  too  many, 
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the  time  step  duration  should  be  shortened. 


The  convergence  criterion  for  the  entire  solving  process  is 


max 


l  u. 


<  «> 


(20) 


In  the  equation  •,  —  io~* 

When  the  initial  field  is  under  consideration  in  the  paper, 
the  physical  values  at  the  boundary  points  are  taken  as  given 
values.  Thus,  the  inherent  boundary  condition  of  step  5  becomes 


Aar  —  0 


(21) 


VI.  Computation  Examples  and  Discussion 

Since  the  convection  terms  in  Eq.  (8)  are  not  separately 
integrated,  mention  of  the  natural  boundary  conditions  of  the 
exit  boundary  is  appropriate.  In  the  following  computation 
examples,  all  conditions  at  the  exit  boundary  are  taken  as 

?--0 ,  (22) 

Ox 

1.  Reflection  of  nonviscous  shock  wave  at  solid  wall 
As  one  of  the  computation  examples,  in  the  paper  the  shock 
wave  reflection  problem  is  the  first  to  be  calculated  (Mach 
number  of  incoming  flow  -  2.9,  and  the  included  angle  between 

O 

the  incident  direction  and  the  horizontal  direction  is  29  ).  In 
the  computations,  the  number  of  elements  is  44  x  11;  the 
homogeneous  flowfield  is  the  initial  condition. 


Fig.  4.  Equal-Mach-number  curves  for  shock 
reflection  at  solid  wall  w.- 1.»,  aw- 0.06M 


Fig.  4  shows  the  distribution  of  equal-Mach-number  curves 
found  by  computation;  the  figure  clearly  exhibits  the  incident 
and  the  reflection  shock  waves.  Fig.  5  shows  the  pressure 
distribution  at  y  =  0.5  given  in  the  computations.  In  three 
zones,  all  the  results  in  the  paper  are  quite  consistent  with  the 
accurate  solution.  The  shock  wave  positions  from  the 
computations  are  slightly  to  the  rearward  compared  to  those  in 
reference  [7]. 
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Fig.  5.  Pressure  distribution  Fig.  6.  Region  of  computation 
at  y  *  0,5  KEY:  *  -  conditions  of  incoming 

0  -  results  given  in  the  paper  flow 

_  -  reference  [8] 


2.  Flat-plate  laminar  flow 
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Fig.  6  shows  the  region  of  solution.  In  the  combining  layer 
zone  at  the  front  fringe  of  the  flat  plate,  the  velocity  drift 
and  temperature  jump  are  neglected.  Thus,  the  boundary 
conditions  of  the  entire  wall  surface  at  y  =  0  can  be  written  as 

*  -  -  o  )  r 


T  —  Tm  —  1  + 


(23) 


In  the  equation,  the  wall  surface  temperature  is  taken  as  the 
total  temperature  of  the  free  incoming  flow. 

Fig.  7  shows  the  pressure  distribution  at  the  wall  surface. 
We  can  see  from  the  figure  that  the  results  of  computation  in  the 
paper  are  relatively  consistent  with  the  results  given  by  the 
weak  interference  theory  in  reference  [8];  in  addition,  there  is 
no  fluctuation  phenomenon  in  the  numerical  solution  at  the  front 
fringe  of  the  flat  plate.  As  revealed  in  the  numerical 
experiments,  at  the  front  fringe  of  the  flat  plate  when  the  step- 
duration  Reynolds  number  Re^x  <  50,  there  is  no  fluctuation  in 
the  numerical  solution. 

Fig.  8  shows  the  velocity  and  temperature  profiles  at 
x/L  *»  1.0.  From  the  figure,  we  can  see  that  the  boundary  layer 
is  approximately  at  y/L  -  0.25,  and  the  front-edge  shock  wave  is 
approximately  sited  at  y/L  «*  0.57.  This  conclusion  is  to  be 
expected.  x/L  «  1.0  is  sited  downstream  of  the  front-fringe 
combined  layer.  Strong  interference  disintegrates  as  weak 
interference.  In  Fig.  8(b),  the  negative  temperature  gradient  is 
exhibited  in  the  neighborhood  of  the  wall  surface.  This 
indicates  heat  absorption  is  needed  from  the  wall  surface.  Given 
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by  Eq.  (23),  the  wall  surface  temperature  is  the  total 
temperature  of  free  flow  as  the  adiabatic  wall  temperature  when 
Pr  =  1.  However,  in  the  case  of  air  with  Pr  **  0.72,  the 
restoration  coefficient  R  is  always  smaller  than  1.  In  other 
words,  the  adiabatic  wall  temperature  should  be  smaller  than  the 
wall  temperature  given  in  Eq.  (23).  Therefore,  in  this  case  the 
flat  plate  should  serve  as  the  adiabatic  wall  surface. 


<»)  “  1000  (b'<  Re.t  -  6000 

Fig.  7.  Pressure  distribution  at  wall  surface 
LEGEND:  0  -  results  in  the  paper 
_  -  reference  [9] 


k/U-  -  tit. 

Fig.  8.  Parameter  profile  for  the  case  when 
Re oo l  “  1000  at  x  ■  0.1 
LEGEND:  o  -  results  in  the  paper 
-  reference  [10] 

(a)  Velocity  distribution  (b)  -  Temperature 
distribution 

KEY:  1  -  shock  wave  2  -  boundary  layer 
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We  also  see  from  Fig.  8  that  the  velocity  distribution  from 
the  computations  in  the  paper  are  in  relatively  close  agreement 
with  the  values  in  [9].  However,  the  temperature  distribution 
shows  certain  discrepancies  within  the  boundary  layer.  On  the 
one  hand,  this  is  because  the  distribution  rule  of  wide-direction 
density  in  the  vicinity  of  the  wall  surface  described  in 
reference  [9]  makes  specifications;  on  the  other  hand,  this  is 
because  only  the  viscosity  terms  that  are  larger  than  the 
magnitude  of  OfRc;'^)  are  retained  in  the  quasi-simplified 
Navier-Stokes  equation  set  used  in  this  paper.  Obtained  from 
this  form  of  the  equation  set,  the  temperature  distribution  in 
the  vicinity  of  the  wall  surface  often  shows  larger  discrepancies 
from  the  solution  obtained  from  the  entire  Navier-Stokes  equation 
set. 

3.  Supersonic  compression  corner  laminar  flow 

The  foregoing  method  is  used  in  the  paper  to  calculate  the 
flow  at  the  two-dimensional  compression  corner  flow  of  »  3 
and  Rc..t  —  6  x  10*  .  Fig.  9  shows  the  region  of  solution  of  the 

problem;  the  boundary  conditions  of  the  wall  surface  are  shown  in 
Eq.  (23). 

In  the  paper,  the  position  of  the  left  boundary  is  located 
at  x  ■  0.3571;  the  physical  value  at  the  boundary  is  taken  from 
the  computation  result  of  flat-plate  flow  mentioned  previously. 
The  number  of  matched  points  is  35  x  26;  the  matching  mesh  along 
the  x  and  y  directions  is  equally  spaced,  respectively,  at 
x  «  0.0428  and  y  «  0.0035. 
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Fig.  9.  Region  of  solution  and  boundary  for 
compression  corner  supersonic  flow 
KEY:  1  -  Front  edge  shock  wave  2  -  Extrapola¬ 
tion  of  simple  wave  3  -  Upstream  condition  is 
based  on  the  flat-plate  solution  in  the  paper 
4  -  Downstream  condition  5  -  Wall  surface 
condition 


* 


Fig.  10.  Wall  surface  pressure  distribution 
when  (O  ■  7.5° 

LEGEND:  0  -  results  in  the  paper 

_  -  r  .erence  [10]  -  -  -  -  nonviscous 

accurate  solution 


Fig.  10  shows  the  wall  surface  pressure  distribution 
obtained  from  computation  when  (0  »  7.5°.  For  comparison,  the 
figure  also  shows  the  numerical  results  of  Carter’s  difference 


16 


method  [9].  From  the  figure  we  can  see  that  the  pressure 
continues  to  drop  within  a  certain  distance  beginning  from  the 
left  boundary.  This  indicates  that  the  influence  of  the  upstream 
has  been  taken  into  account.  Because  of  the  loss  due  to  a  shock 
wave  of  a  nonviscous  flow,  the  wall  surface  pressure  of  the 
downstream  region  thus  calculated  is  higher  than  the  nonviscous 
flow  pressure. 

VII.  Conclusions 

In  the  paper,  the  SUPG  finite-element  method  is  developed 
and  applied  to  computations  of  supersonic  viscous  flow  in 
obtaining  relatively  satisfactory  numerical  results,.  From  the 
computations  in  one-dimensional  and  two-dimensional  problems,  the 
SUPG  method  certainly  can  restrict  the  fluctuations  in  a  shock 
wave  both  downstream  and  upstream*  However,  there  is  a  problem 
of  selecting  the  free  parameters.  Therefore,  it  is  necessary 
that  we  further  develop  the  finite-element  format  of  nonfree 
parameters  that  do  not  show  fluctuations. 

The  first  draft  of  the  paper  was  publicly  read  at  the  Fifth 
Session  of  the  All-China  Computational  Fluid  Mechanics  Conference 
in  April  1990,  at  Huangshan  Mountain.  Xu  Guoqun  is  currently 
working  at  the  Jiangsu  Provincial  Investment  Corporation.  The 
paper  was  received  for  publication  on  4  September  1990. 
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